Generación de números aleatorios
Marcos García González (h[e]rtz)
Date: Verano 2004
Documento facilitado por la realización de la asignatura Métodos informáticos de la física de segundo curso en la universidad autónoma de Barcelona en el transcurso del año 2003-2004.
Un problema básico que nos encontramos habitualmente es el de obtener secuencias de números uniformemente distribuidos en un intervalo
.
La diferentes posibilidades para resolver dicho problema son:
i) Buscar en tablas de números aleatorios publicadas (libros, internet ...);
ii) Observar un proceso físico tal como la desintegración radiactiva, el ruido eléctrico ...;
iii) Los lenguajes de programación y las hojas electrónicas incluyen una función para generarlos
iv) Mediante algorismos de generación de números aleatorios
Las principales ventajas de los generadores de números aleatorios son:
- Rapidez
- Comodidad
- Reproducibilidad
- Portabilidad
Y la desventaja fundamental:
- Las secuencias obtenidas no son realmente aleatorias, ya que se obtienen con operaciones deterministas. Solo podemos obtener secuencias pseudo-aleatorias, que a su vez satisfacen algunos criterios de aleatoriedad adecuados.
Los números generados deben cumplir ciertas características para que sean válidos. Dichas características son:
1. Uniformemente distribuidos.
2. Estadísticamente independientes.
3. Su media debe ser estadísticamente igual a 1/2.
4. Su varianza debe ser estadísticamente igual a 1/12.
5. Su periodo o ciclo de vida debe ser largo.
6. Deben ser generados a través de un método rápido.
7. Generados a través de un método que no requiera mucha capacidad de almacenamiento de la computadora.
Normalmente se utilizan números enteros, ya que su aritmética es exacta y rápida. Se generan enteros
entre 0
y
, y
da valores reales en el intervarlo
.
En general los algoritmos utilizan relaciones de recurrencia del tipo
en el caso de recurrencia simple, o bien
para el caso de una recurrencia de orden
![$ k$](img8.png)
.
Se necesitará dar un
valor inicial para comenzar el algoritmo (
![$ k$](img8.png)
valores para recurrencias de orden
![$ k$](img8.png)
).
Estos generadores son los más utilizados y los más conocidos. Se basan en la relación de recurrencia
donde
![$ a$](img10.png)
es el multiplicador y
![$ m$](img11.png)
el módulo.
- Hay
![$ m$](img11.png)
valores posibles de
![$ N_i$](img2.png)
, entre 0 i
![$ m-1$](img12.png)
.
- La secuencia es
periodica: cuando vuelve a aparecer un número por segunda vez, la secuencia se vuelve a repetir. El periodo depende de los valores de
![$ a$](img10.png)
,
![$ c$](img13.png)
y
![$ m$](img11.png)
, así como del valor inicial; nótese que el máximo posible es
![$ m$](img11.png)
.
Recordemos que lo que nos interesa para trabajar con un buen generador de números aleatorios es que la distribución de los números obtenidos tiene que ser uniforme, no deben de haber correlaciones entre los terminos de la secuencia, el periodo debe ser lo más largo posible, y el algorismo debe ser de ejecución rápida.
Las limitaciones más importantes de los generadores son su periocidad (normalmente el periodo no suele ser más grande de
![$ 2^{32} \approx 4 \times 10^9)$](img14.png)
y la posible presencia de correlaciones entre términos consecutivos de la secuencia. Una manera sencilla de suprimir éstas limitaciones es
![$ ''$](img15.png)
desordenar
![$ ''$](img15.png)
un poco la secuencia mediante el siguiente procedimiento:
Se parte de un generador que da enteros aleatorios entre
0
y
![$ m-1$](img12.png)
, y en primer lugar se genera con el GCL un vector que contiene una lista de
![$ N$](img16.png)
enteros aleatorios
![$ j_n$](img17.png)
, así como un entero aleatorio
![$ y$](img18.png)
. Se determina el índice
![$ k = [y * N/m]$](img19.png)
, entre
0
y
![$ N-1$](img20.png)
.
El elemento
![$ j_k$](img21.png)
de la lista se da como un nuevo nombre aleatorio, y se reasigna a la variable
![$ y$](img18.png)
el valor
![$ j_k$](img21.png)
. El valor de
![$ j_k$](img21.png)
se renueva con el GCL, y se vuelve a repetir los pasos desde la determinación del índice
![$ k$](img8.png)
.
En estos generadores cada nuevo número entero aleatorio
![$ N_i$](img2.png)
, se obtiene manipulando los bits del número anterior,
![$ N_{i-1}$](img22.png)
. En lenguaje C, esto se puede hacer facilmente utilizando operadores sobre bits,
![$ », «, \wedge, \vert, \&$](img23.png)
.
Las grandes ventajas de estos generadores es que son generadores muy rápidos que tienen un periodo muy largo. La fomentación teorica en la que se basan es diferente a la de los GCL. Los generadores de Fibonacci se basan en una recurrencia del tipo
donde
![$ r<s$](img25.png)
son enteros dados y
![$ \circ$](img26.png)
denota alguna de las operaciones
![$ +, -, \times, \wedge$](img27.png)
.
Este tipo de generador precisa iniciar (con otro generador) y mantener una lista de los últimos
![$ s$](img28.png)
números generados.
Otros tipos de generadores los podemos encontrar en:
- W.H. Press, S.A. Teukolski, W.T. Vetterling i B.P. Flannery,
Numerical Recipes in C, Cambridge University Press.
-D.E. Knuth,
The Art of computing programming, 2: Seminumerical Algorithms, Addison-Wesley.
Antes de aceptar un nuevo generador hace falta verificar que satisface una série de pruebas, lo que llamaremos
pruebas de aleatoriedad.
Éstas consisten basicamente en realizar dos tipos de pruebas, empíricas y teoricas. Si los generadores superan estas pruebas, podremos asegurar que estamos ante un generador de números aleatorios bastante competente.
Pruebas empíricas (sobre la muestra de la secuencia)
- Test de uniformidad: hace falta que los valores esten uniformemente distribuidos en
![$ [0,1]$](img1.png)
. Se puede realizar un test
![$ \chi^2$](img29.png)
. Alternativamente, podemos estimar los momentos de orden
![$ k, X^k=1/N \sum_i (x_i)^k$](img30.png)
y comprobar que se aproximan a sus correspondientes valores teoricos,
![$ 1/(k+1)$](img31.png)
.
- Test serial: se generan parejas de valores
![$ (x_1,x_2)$](img32.png)
, y se comprueba si se distribuyen uniformemente en el cuadrado
![$ [0,1]\times [0,1]$](img33.png)
.
- Test de correlaciones: se determina la correlación entre números separados
![$ k$](img8.png)
lugares en la secuencia,
![$ C(k) = 1/N \sum_i x_i x_{i+k}$](img34.png)
. Su valor tendría que acercarse a 1/4.
Pruebas teoricas (sobre toda la secuencia): La sencillez de los generadores de congruencia lineal permiten demostrar propiedades importantes:
- Para determinar valores de
![$ a, c$](img35.png)
i
![$ m$](img11.png)
se obtienen secuencias de periodo máximo
![$ m$](img11.png)
. Por ejemplo, si
![$ m$](img11.png)
es una potencia de 2, bastará con que
![$ c$](img13.png)
sea impar y
![$ a$](img10.png)
sea igual a un múltiple de 4 mas 1.
- Test espectral: si se forman vectores con
![$ k$](img8.png)
valores consecutivos,
estos forman hiperplanos paralelos en el espacio k-dimensional. La separación
![$ d_k$](img37.png)
entre los planos tiene que ser la mínima posible.
El problema a tratar será el de obtener una secuencia con densidad de probabilidad dada
![$ w(y)$](img38.png)
, definida en el intervalo
![$ (y_{min}, y_{max})$](img39.png)
a partir de una secuencia de números con distribución uniforme
![$ U(0,1)$](img40.png)
.
Para resolver este problema utilizaremos el
método del cambio de variable. Que concretamente consistirá en buscar una transformación
![$ y = f(x)$](img41.png)
para obtener la distribución deseada. Si la densidad de probabilidad de
![$ x$](img42.png)
es
![$ p(x)$](img43.png)
, tenemos:
Si
![$ x \in U(0,1)$](img46.png)
la integración de esta ecuación resulta ser
Si sabemos calcular la integral, obtenemos una relación del tipo
![$ X = g(Y)$](img48.png)
, que se tiene que invertir para poder obtener
Distribución exponencial.
Para obtener
![$ Y$](img50.png)
con distribución
![$\displaystyle w(y) = \left\{\begin{array}{ll} 0, & y<0
\lambda e^{-\lambda y}, & y \in [0, +\infty)
\end{array} \right.$](img51.png) |
(1) |
entonces nos queda la ecuación
por lo tanto, el cambio adecuado será
Distribución normal. Para obtener
![$ Y \in N(0,1)$](img54.png)
, nos queda la siguiente ecuación
La integral no se puede resolver analíticamente, y menos aun invertir la relación para obtener
![$ Y = f(X)$](img56.png)
.
La alternativa es aplicar el
algoritmo de Box-Muller. Que permite obtener 2 valores independientes
![$ x_1$](img57.png)
,
![$ x_2$](img58.png)
con distribución
![$ N(0,1):$](img59.png)
si
![$ x_1, x_2$](img60.png)
son las coordenadas de un punto del plano, la densidad de probabilidad de sus coordenadas polares es
Por lo tanto,
![$ \theta \in U (0,2\pi)$](img62.png)
i
![$ r$](img63.png)
tiene la densidad de probabilidad
![$ re^{-r^2/2}$](img64.png)
. Se pueden obtener a partir de 2 números
![$ z_1, z_2 \in U(0,1)$](img65.png)
con las transformaciones
![$ \theta = 2 \pi z_1, r = \sqrt{-2 \textrm{ln}z_2}$](img66.png)
. Por lo tanto
![$ x_1, x_2$](img60.png)
se pueden obtener con las transformaciones
![$\displaystyle \left\{\begin{array}{ll} x_1 = \sqrt{-2 \textrm{ln}(z_2)} \te...
...
x_2 = \sqrt{-2 \textrm{ln}(z_2)} \textrm{sin}(z_1)
\end{array} \right.$](img67.png) |
(2) |
Método de aceptación-rechazo: Es un método sencillo y general, aunque en ocasiones no muy eficiente.
Sea
![$ p(x)$](img43.png)
definida en un intervalo finito,
![$ (a,b)$](img68.png)
, y
![$ M$](img69.png)
una cota superior de
![$ p(x)$](img43.png)
. Se generan dos valores aleatorios
![$ x \in U (a,b)$](img70.png)
y
![$ p \in U(0,M).$](img71.png)
Entonces:
![$\displaystyle \left\{\begin{array}{ll} \textrm{se acepta} x \textrm{si} ...
... p
\textrm{se rechaza} x \textrm{si} p(x) < p
\end{array} \right.$](img72.png) |
(3) |
De esta manera el valor
![$ x$](img42.png)
aparece con una densidad de probabilidad
![$ p(x)$](img43.png)
, aunque no se aprovechan todos los valores que obtenemos en la realización del método. De todos modos, este es el único método disponible cuando la distribución de probabilidades es complicada. De hecho es la base del
algoritmo de Metropolis, posiblemente el más utilizado en física computacional.
DoFinal(); ?>